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Abstract 

Beta regression models provide an adequate approach for modeling continuous 
outcomes limited to the interval (0, 1). This paper deals with an extension of beta 
regression models that allows for explanatory variables to be measured with er- 
ror. We follow the structural approach, in which the covariates measured with 
error are assumed to be random variables. We present three estimation methods, 
namely maximum likelihood, maximum pseudo-likelihood and regression cali- 
bration. Monte Carlo simulations are used to evaluate the performance of the 
proposed estimators and the naive estimator. Also, a residual analysis for beta re- 
gression models with measurement errors is proposed. Our results are illustrated 
in a real data set. 

Keywords: Beta regression model; Errors-in-variables model; Gauss-Hermite 
quadrature; Maximum likelihood; Maximum pseudo-likelihood; Regression cal- 
ibration. 



1. Introduction 

Errors-in-variables models, also called measurement error models, are widely 
applicable in many research areas since they allow for the presence of explicative 
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variables that are measured with errors or that cannot be observed directly (latent 
variables). Many examples and applications of these models are considered in the 
books by |Fuller| ( [T9871 ), [Carroll et al. \ ( |2006l ) and [Cheng and Van Ness| ( fT999l ). It 
is well known that measurement errors cause biased and inconsistent parameter 
estimates and lead to erroneous conclusions in inferential analysis. Errors-in- 
variables models are specified in such a way that the distribution of the response 
variable, y, is assumed to depend on covariates, x, which are imprecisely mea- 
sured, and observable variables, w, are seen as surrogates for the unobservable 
true covariates. The classical linear errors-in-variables model has been extensively 
discussed in the literature, particularly under the normality assumption for the dis- 
tribution of the unobservable variables. For a systematic review of such models 
see 



Fuller] ([1987 ) and [Cheng and Van Ness[ ( |T999| ); see also [Arellano- Valle and 



Bolfarine[ ( [T996i ) and [Castro etaL\^200E^ . 

The beta regression models provide an adequate approach for modeling con- 
tinuous outcomes limited to the interval (0, 1), or more generally, limited to any 



open interval (a, 6) as long as the limits are known (Ferrari and Cribari-Neto 



2004). Although the literature on beta regression has grown fast in the last few 
years, errors-in-variables models with beta distributed outcomes is an unexplored 
area. 

A beta regression model assumes that the response variable, y, has a beta 
distribution with probability density function 



m 



r(/x0)r[(i-/x)0] 



y 



{i-m)</'-i 



0<y<l, 



(1) 



where r(-) is the gamma function, < < 1 and > 0, and we write y ~ 
Beta(/i, 0). Here, fi = E(?/) and is regarded as a precision parameter since 
Var(?/) = — + 0). For independent observations yi,y2, . ■ . , yn, where 

each yt follows a beta density Q with mean /it and unknown precision parameter 
0, the beta regression model defined by Ferrari and Cribari-Neto (2004) assumes 
that 



(2) 



with a. G W"' being a column vector of unknown parameters, and with zj = 
{zti, . . . , Ztp^) being a vector of pa fixed covariates (pa < n). The link function 
g{-) : (0, 1) — 7- M is assumed to be a continuous, strictly monotone and twice 
differentiable function. There are many possible choices for g{-), for instance, the 



2 



logit link, = log[/it/(l - yUt)], the probit link, g{ij,t) = ^ ^{^t), where $(■) 
is the cumulative distribution function of the standard normal distribution, and the 
complementary log-log link, g{iJit) = log[— log(l — fit)]. 



Extensions for the beta regression model proposed by Ferrari and Cribari 



Neto| ( |2004| ) that allow the precision parameter to vary across observations, or that 
involve non-linear structures for the regression specification of the mean and the 
precision parameter, are presented by jSmithson and Verkuilen| ( [2006| ), |Simas et al. 



( |2010 ), among others. The beta regression model with linear specification for the 
transformed mean and precision parameter is given by ([T]), ([2]) and 

h{<Pt) = v77, (3) 

where 7 G (pa + < n) is a column vector of unknown parameters, = 
{vti, • ■ ■ , Vtp^y is a vector of fixed covariates, h{-) : (0, 00) — > M is a strictly 
monotone, twice differentiable link function. A possible choice for h(-) is = 

log(0t)- 

The purpose of this paper is to extend the beta regression model ([l])-([3]) to the 
situation where some covariates are not directly measured or are measured with 
error. A practical application of errors-in-variables beta regression models will be 
illustrated in a study of the risk of coronary heart disease as a function of low- 
density lipoprotein (LDL) cholesterol level ("bad cholesterol") and body mass 
index (BMI). The dataset consists of observations of systolic blood pressure 
(SBP), diastolic blood pressure {DBF), BMI and total cholesterol level (TC) 
in a group of 182 smoking women aged 50 to 87 years. The total cholesterol 
may be considered as a surrogate of LDL, which is a covariate of interest, and 
whose direct measure is more expensive and time consuming. The difference 
between SBP and DBF results in what is known as the pulse pressure, PP = 
SBP-DBP, and the relative pulse pressure is RPP = {SBP-DBP)/ SBP = 
PP/SBP. Small values of RPP, RPP < 25% say, is indicative of risk of heart 
disease ( [American College of Surgeons} |2008[ p. 58). Notice that the response 



variable, RPP, is continuous and limited to the unit interval, and that one of the 
covariates, namely LDL, is not measured directly. 

This paper is organized as follows. In Section [2[ we present an errors-in- 
variables beta regression model under the structural approach, and the correspond- 
ing likelihood function. In Section |3j we present three different estimation meth- 
ods, namely maximum likelihood, maximum pseudo-likelihood, and regression 
calibration. In Section]?} we perform a simulation study to evaluate and compare 
the performance of the three estimation approaches. In Section ]5} we propose a 
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residual analysis. Section [6] presents an application of the proposed model. Con- 
cluding remarks are presented in Section |7] 

2. Model and likelihood 

Let yi, . . . ,yn be independent observable random variables arising from a 
sample of size n, such that yt has a beta distribution with probability density func- 
tion ([T]) with parameters = ^it and (j) = (j)t. In the following, we assume that Ht 
and 0t may depend on covariates and unknown parameters. In practice, some co- 
variates may not be precisely observed, but, instead, may be obtained with error. 
The model considered in this paper assumes a linear structure for the specification 
of the mean and the precision parameters, and also assumes that both specifica- 
tions may involve covariates measured with error. Specifically, we replace the 
mean submodel ([2]) and the precision submodel ([3]) by 

g{yit) = zja + ^jf3, (4) 
hi(l>t) = v77 + m7A, (5) 

respectively, where f3 E Wf^ , A G are column vectors of unknown parame- 
ters, Xi = {xti, ■■■ , xtp^y and mt = {rriti, ■ ■■ , mpj~^ (Pa + Pfi + + P\ < n) 
are unobservable (latent) covariates, in the sense that they are observed with error. 
The vectors of covariates measured without error, and v^, may contain variables 
in common, and likewise, and nif Let Sj be the vector containing all the unob- 
servable covariates. For t = 1, . . . ,n, the random vector is observed in place 
of St, and it is assumed that 

wt = To + Ti o St + et, (6) 

where is a vector of random errors, tq and ti are (possibly unknown) parame- 
ter vectors and o represents the Hadamard (elementwise) product. The parameter 
vectors tq and ri can be interpreted as the additive and multiplicative biases of 
the measurement error mechanism, respectively. If Tq is a vector of zeros and ri 
is a vector of ones, we have the classical additive model = Sj + e^. Here, we 
follow the structural approach, in which the unobservable covariates are regarded 
as random variables, i.e. we assume that Si, . . . , s„ are independent and identi- 
cally distributed random vectors. In this case, it is also usual to assume that they 
are independent of the measurement errors ei, . . . , e„. Moreover, the normality 
assumption for the joint distribution of and is assumed. The parameters of 
the joint distribution of and St is denoted by d. 
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Let (i/i, wi), . . . , {nn, w„) be the observable variables. We omit the observable 
vectors zt and in the notation as they are non-random and known. The joint 
density function of (yt,wt), which is the observation for the t-th individual, is 
obtained by integrating the joint density of the complete data (yt, w^, s^), 

fiVu wt, Si; e, 6) = /(yt|wt, St; e)f{st, wf, 6), 

with respect to s^. Here, 6 = (ck^, /3^, 7^, A^)^ represents the parameter of 
interest, and 6 is the nuisance parameter. The joint density /(wj, St;S), which 
is associated to the measurement error model, can be written as /(w^, St;S) = 
f(wt\st; S)f{st\S) as well as f(wt, s*; S) = /(st|wt; d)f{wt\d). In this work we 
assume that, given the true (unobservable) covariates St, the response variable yt 
does not depend on the surrogate covariates Wt; i.e. f(yt\wt, s^; 0) = f(yt\st; 6). 
In other words, conditionally on s^, y^ and Wj are assumed to be independent 



(Bolfarine and Arellano- Valle 1998 1. Therefore, the density function of (y^, w^) 



is given by 

/oo /"OO 
••• / f{yt,^t,St;0,S)dst, 
-00 J —00 

00 poo 



f{yt\st;0)f{wt,St;5)dst. 
The log-likelihood function for a sample of n observations is given by 

^ /"OO /"CxD 

eie,S) = Y,^og ••• / f{yt\sue)f{st\wt;S)f{wt;6)dst, 

J —00 J —00 

^ poo /"OO 

= ^^ogf{wt;S) + ^log •••/ f{yt\st;e)f{st\wt]S)dst. 
,_i .-I J— 00 J —00 



(7) 



t=l i=l 



In general, the likelihood function involves analytically intractable integrals and, 
hence, approximate inference methods need to be considered. In the next section, 
we present three different approaches to estimate the parameters. 

In order to facilitate the description of the estimation methods, we assume that 
a single covariate, Xt, is measured with error, and that it is used for the specifica- 
tions of both the mean and precision submodels. We then have pj3 = p\ = 1 and 
St = Xt = TUt. We also assume independence and normality of random errors. 
The methodologies presented in this paper can be extended to the situation where 
and are distinct, or when covariates measured with error appear only in the 
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specification of the mean or the precision parameter. To be specific, from now on, 
the model under consideration is summarized as follows. For i = 1, . . . ,n. 



yt\xt, wt ~ Beta(yUt, (pt), (8) 
gifJ-t) = zja + xtl3, h{(f)t) = vjj + XtX, (9) 

Wt = To + nxt + et, xt '~ A^(/ix, cr^), et '~ A^(0, a^), (10) 

with Xt and Cf, for t,t' = 1, . . . ,n, being independent. The unknown parameter 
vectors a. and 7 were defined above, and /3 G M, A G M, /i^: G M and cr^ > 
are unknown parameters. Note that it is assumed that the conditional distribution 



of yt given {xt,Wt) does not depend on Wf. Also, if tq = and ti = 1, (10) 
corresponds to the classical additive error model Wt = Xt + Ct. 



From ( 10) we have 



ind -^j, 2 2 2\ i ind 2\ 

~ N(ro + ri/i^,ri(T^, + (Tj, ~ N(/i^,|„,, a^^|^J, (11) 



where 



/J-Xtlmt = l^x + K[Wt - (to + TijJ,^)], (T^^i^^ = (Tg/c^/ri, (12) 

with kx = ricr^/(rfcr^ + cTg) being known as the reliability ratio. To avoid non- 
identifiability of parameters we assume that (ro,ri,crg) or {To,Ti,kx) are either 
known parameters or they are estimated from supplementary information, typi- 
cally replicate measurements or partial observation of the error-free covariate. In 
any case, either of these vectors are regarded as known quantities in the inferential 
procedure. Hence, the nuisance parameter vector is <5 = (/i^., cr^)^. 

The log-likelihood function in (|7]) for n observations taken from the model 
described in ([8]), ([9]) and ( 10 ) is given by 



(0,5) = ^£i,((5) + ^£2t(^,<5), (13) 



t=i t=i 
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where 



-log[2vr(r,a, + aj] - 2(rM + a,^) ^ 



oo 



oo 



1 



exp 



xt\wt 



(14) 

cia;t.(15) 



3. Estimation 



3.1. Maximum likelihood estimation 



The second term of the log-likelihood function £(0, 5) in (13 1, which depends 



on a non-analytical integral (as can be seen in ( 15 )), can be approximated using 



the Gauss-Hermite quadrature, which consists of the approximation 



/oo Q 



(16) 



where rjg and Ug represent the g-th zero and weight, respectively, of the orthogonal 
Hermite polynomial of order Q (number of quadrature points); see, for instance. 



(xt-/i.,I^J/j2a2^l^^ in( 



Abramowitz and Stegun| ( [r972[ Section 20.4.46). Using the transformation Uf 

15 1, we have that xt = fixt\wt + J'^'^lA^.'^t and dxt 



2a^^^^^dut. Hence, by applying (16) in (15), the log-likelihood function (|13|) 
can be approximated by 

n n / Q 



t=i 



t=i 



.9=1 



where 



= iogr(0)-iogr(/i0)-iogr[(i-/i)0] 

+ - 1) \og{yt) + [(1 - /i)</. - 1] log(l - yt), 



(17) 



tq) 



zja + x;i3, h{(f)tq) =vJ-f + x;X, xl = + 2al^^^r]q 



where and are given in \\2\. 
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For n and Q sufficiently large, the distribution of (0 ,6 Y , the estimator 
of {6^ ,5^Y obtained by solving the system of equations d£a{0,5)/d9 = 0, 
dia{6, d)/d6 = is approximately normally distributed with mean {6'^,6~^Y 
and covariance matrix J^^iO, S), where 



Ja(6',5) 



(Guolo, 2011 ). For computational implementation, the derivatives of ia{6, 6) with 



respect to the parameters can be analytically obtained or numerical derivatives can 
be used. Our numerical results were obtained using numerical derivatives. 

3.2. Maximum pseudo-likelihood estimation 

The central idea of the maximum pseudo-likelihood estimation method is to 
replace the nuisance parameters by consistent estimates in the log-likelihood func- 
tion ( [T3| ). The resulting function can be regarded as a pseudo-log-likelihood func- 
tion that depends on the parameters of interest only ( Guolo[ 2011 ). 

The log-likelihood function ( [T3] ) is maximized in two steps. First, we esti- 
mate the nuisance parameter vector <5 by maximizing the reduced log-likelihood 
function 



\5) = Y,^n{5), 



t=x 



(18) 



where ^lt{S) is given in ( 14). Second, the estimate 8 obtained from the maximiza- 
tion of ( [T8] ) is inserted in the original log-likelihood function (13), which results 
in the pseudo-log-likelihood function 



(19) 



t=i 



i=l 



As in i{0, S), the second term in £p{0; S) is analytically intractable. However, 
unlike the integral in i{0 , S), the integral in £p{0; 6) depends on the parameter of 
interest only. From (16), it is possible to approximate i2t{d, <5) by a summation. 



After some algebra, we have that an approximate pseudo-log-likelihood for the 



beta regression model with one covariate measured with error is given by (19) 



8 



with i2t{0,S) replaced by 



-2t 




exp 



liJitg, (f>tq) 



where 0) is given in ( 17 1, 



9{%q) = ^t " + K(t^tq) = Vi 7 + XtA, xt = + ^2a\^^r]q, 

with and being estimates of yU^-^i^j and respectively, obtained 

from the maximization of ( [T8| ). 

It can be shown that, under regularity conditions ( Gong and Samaniego[ 1981 1, 
the asymptotic distribution of ^Jn{0 — 0) is normal with mean zero and covariance 
matrix 



I 



ee 



66' 



where 
^66 



t=l 



dXt{6,6) 
0686^ 



t=i 

n 

E 

t=i 



86 



86 



8irtiS) f8US) 



8S 



86 



Ee 



t=l 



'66 



8Hpti6,6) 



8686^ 



E 



dU6) f 8ipt{6,6) 
86 V 86 



with 6 replaced by 6, and with irt{6) and ipt{6, 6) being the t-th element of the 
log-hkelihood functions £ri6) and ip{6, 6) in ([Tsj) and ([T9|), respectively (|Guolo 



201 1| ). Further details can be found in [Carroll et p006| , Section A. 6. 6). For 



the errors-in- variables beta regression model considered here, Iq0, Igg and I^^ 
do not have closed form. We suggest to replace the expected information ma- 
trix by the observed information matrix. For computational implementation, the 
needed derivatives can be analytically or numerically obtained. We used numeri- 
cal derivatives in our applications. 



9 



3.3. Regression calibration estimation 

The regression calibration method has been widely used in errors-in- variables 
modeling due to its simplicity; see [Carroll et a/.| ( [2006[ Chap. 4), [Freedman et al. 
( |2008| ), [TTiurston et ( [2005] ) and |Guolo| ( |20TT] ). The central idea is to replace 
the unobservable variable, xt, by an estimate of the conditional expected value of 
Xt given wt, 'E{xt\wt), in the likelihood function. Let r(wt, 6) = 'E{xt\wt) be the 
calibration function. The replacement of the unobservable covariate Xt by r(wt, 6) 
establishes a modified model for the data. Here, 6 is an estimate of 6. 

For our errors-in-variables beta regression model, the calibration function is 
r{wt,6) = ^ixt\wt as defined in (12) 
EILi wt/n and si = ^"^i 



[Wt 



w 



From (|11|) and (|12|) 

1 



we have that w = 
are optimal estimates of tq + ri/i^ 



and rf (T^+cXg , respectively. These estimates can be used to estimate the calibration 
function. 

By inserting the calibration function in the conditional density function of yt 
given Xt, we obtain the modified log-likelihood function 



(6) 



t, Vt) 



(20) 



t=i 



where /(/i, 0) is given in ( 17 ), 



g{lj,t) = zj(x + xt(3, h{(j)t) = J + XtX, 



with Xt being the estimated calibration function. Note that the modified log- 
likelihood function in ([20]) only depends on the parameter of interest, 6. The 



regression calibration estimate of is obtained from the system of equations 
d£rc{G)/dO = 0, which requires a numerical algorithm. Since irc{G) coincides 
with the log-likelihood function for the usual beta regression model, Xt acting 
as an observable covariate, these equations can be numerically solved from avail- 



able computational packages, for instance the betareg package ( iCribari-Neto and 
Zeiles[ 2010) implemented in the R platform. 



We note that (20) is a pseudo-likelihood as presented in Gong and Samaniego 



(1981). Gong and Samaniego ( 1981[ ) present the maximum pseudo-likelihood 
method as follows. Consider a sample of size n with joint density function /(■; 6, 6), 
and assume that an estimate 5 of 5 is available. Hence, the maximum pseudo- 
likelihood estimate, 0, is obtained by maximizing f{-;6,6). One expects to 
have good properties for large samples whenever <5 has; for further details see 
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Delgado|([T995|). From|Gong and Samaniegol ([T98T]), [PEiel (fT986l) and|Delgado 



( 1995 ) it can be shown that, under regularity conditions, the asymptotic distribu- 
tion of ^/n{0 — 6) is normal with mean zero and covariance matrix 



T* — 1 / T* I T* X~l* T*T \ T* — 1 

^ - ^ee vee + ^66^55^65 ) ^ee ' 



where 



^66 



EE 



t=i 

n 



EE 



i=l 



de V de 

'dirct{0,6) fdirct{0,6) 



86 



dS 



t=l 



dS 



dS 



with 5 replaced by S and with irct{6, <5) being the t-th term of irc{6, 6). To obtain 
S* one needs to obtain dircti^^ 6)/ 86, Tqq and I^^. We have 



dirct{0,S) 

dS 
where 



/3{i - nK)Uyl - K)9'{f^t)-' + A(i - nK)ath'{<Pty' 

- Tik^)(pt{y*t - fi;)g'{fit)~^(T-^ + A(l - nk^)ath'{(pty^a-'^ 



«i = l^tiyt - /if) + log(l - yt) - ^((1 - /it)0t) + ipi^t), 



y*t = log 



yt 



i-yt 



(21) 



with ip{-) being the digamma function, i.e., tp{z) = dlogr(2;)/dz for z > 0. 
Here, we are able to obtain Tqq and I^^, since the log-likelihood function does 
not involve integrals. After some algebra, we have 



Z^$UZ ZT<l>Ux Z^CGHV Z^CGHx 
x^^Ux x^CGHV x^CGHx 
VDV V^Dx 
x^Dx 
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where 



Z = (zi, . . . ,z„), V = (vi, . . . , v„), U = diag{Mi, . . . ,M„}, 

$ = diag{0i, . . . , 0J, H = diag{/i'(0i)-\ . . . , /i'(0„)-^}, 
G = diag{^'(/ii)~\ . . . , g'ifin)'^}, C = diag{ci, . . . , c„}, 
D = diagjcii, . . . ,cin}, 



with 



ut = 4>t 



Analogously, we have 



^05 



ZT(U* + HGC*)1„ ZT(Ut + HGCt)l„ 

x^(U* + HGC*)1„ xT(Ut + HGCt)l„ 

VT(D* + HGC*)1„ VT(Dt + HGCt)l„ 

xT(D* + HGC*)1„ xT(Dt + HGCt)l„ 



where 



U*=diagK,...,<}, Ut = diag{nl,...,4}, C* = diag{ct, 
Ct = diag{cl, . . . , 4}, D* = diag{cit, . . . , <}, Dt = diag{4, . , 



with 



4 = c*A;^a^^ rf* = A(l - Tik^)dt, d\ = d^t^a''^, 
and In is a vector of ones(n x 1). 

4. Monte Carlo simulation results 

We now present Monte Carlo simulation results on the performance of the 
different estimation methods described in Section [3l All simulation results are 
based on 5,000 Monte Carlo replications. We consider errors-in-variables beta 
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regression models with log(yUt/(l — fit)) = a + f3xt, and log(0t) = 7 (constant 
precision model) and log(0t) = 7 + Axt (varying precision model), with wt and 
Xt being generated from (10). We set a =2.0, (3 = —0.6, A = 0.5, fXx = 2.5, 
al = 2.7, and 7 = 2.5 for the constant precision model and 7 = 4 for the 
varying precision model. The parameters of the measurement error mechanism 
are assumed to be known, and we set tq = 0, ri = 1, and the following values 
for (Tgi 0.05 and 0.50 (kx = 0.98, and 0.84, respectively). The sample sizes are 
n = 25, 50, 100, 200, and 300. For each simulated sample the parameters were 
estimated under two different settings. First, we ignored the measurement error 
in Xt, i.e. we used what is called the naive method (inaive)', second, we recog- 
nized that Xt is measured with error and estimated the parameters using the three 
methods proposed in this paper: approximate maximum likelihood {£a), (approx- 
imate) maximum pseudo-likelihood (ip), and regression calibration (i'rc)- When- 
ever Gaussian quadrature was required, use set the number of quadrature points 
at Q = 50. The maximization of the relevant (approximate/pseudo/modified) 
likelihoods was performed using the quasi-Newton BFGS nonlinear optimization 
algorithm with numerical derivatives implemented in the function MAXBFGS in 
the matrix language programming Ox ( Doomik[ 201 1 ). 

Figures [T]j4]present plots of the bias and the root-mean-square error of the es- 
timators against sample size, for = 0.98 and 0.84 under the constant precision 
model and the varying precision model. As expected, the naive estimator is biased 
and its bias and mean-square error do not converge to zero as n grows even when 
the reliability ratio is large (i.e. the variance of the measurement error is small 
compared to the variance of x). In other words, the plots suggest that the naive 
estimator is not consistent. For all the cases, the approximate maximum likeli- 
hood and maximum pseudo-likelihood estimators perform similarly. In general, 
their performance is clearly better than that of the regression calibration and naive 
estimators. Under constant precision (Figures [T]j2]), the regression calibration es- 
timator is as biased as the naive estimator for estimating the precision parameter. 
However, for estimating [3, the coefficient associated to the covariate measured 
with error, it performs well if the measurement error variance is small. In fact, 
the regression calibration, approximate maximum likelihood ratio and maximum 
pseudo-likelihood ratio estimators are virtually unbiased for estimating (3 when 
kx = 0.98; also, their mean-square errors converge to zero as n grows. Under the 
varying precision model (Figures |3]j4]), similar conclusions are reached. 

We now turn to the evaluation of confidence intervals constructed from the 
different estimators coupled with their respective asymptotic distributions. Fig- 
ures |5]j6]present plots of the estimated true coverages of confidence intervals con- 
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structed with 95% nominal confidence level, for = 0.98 and 0.84 under the 
constant precision and varying precision models for n ranging from 25 to 300. 
For all the cases, the estimated true coverages of the confidence intervals based on 
the naive estimator decrease as n grows, and hence they cannot be recommended. 
For all the cases, the confidence intervals constructed from the approximate maxi- 
mum likelihood and maximum pseudo-likelihood estimators present true coverage 
close to 95%, more so if the sample size is large. For the constant precision model 
(Figure|5]), the regression calibration estimator produces reliable confidence inter- 
vals for parameters of the mean submodel when the measurement error variance 
is small. However, for estimating the precision parameter, the regression calibra- 
tion estimator produces confidence intervals with true coverage much smaller than 
95% when n is large. Under the varying precision model (Figure |6]), we arrive at 
similar conclusions. 

Overall, we conclude that ignoring the measurement error produces mislead- 
ing inference. Also, if the measurement error variance is small, the regression 
calibration approach is reliable for estimating the parameters of the mean sub- 
model. Moreover, inference based on the approximate likelihood and the pseudo- 
likelihood methods present good performance for the estimation of all the param- 
eters. Since the pseudo-likelihood approach is computationally less demanding 
than the approximate maximum likelihood approach, we recommend the maxi- 
mum pseudo-likelihood estimation for practical applications. 



5. Residual analysis 



Espinheira et al.\ ( |2008[ ) and [Ferrari et al.\ ( |2011| ) proposed the use of stan- 



dardized weighted residuals as a diagnostic tool for beta regression models with 
constant and non-constant precision parameter, respectively. Here, we modify 



the residuals defined by Ferrari et al. (2011 1 to allow for measurement error in 
covariates. We therefore define 



n 



Vt - fJ^t 



vta 



for t = 1, . . . , n, where yl and are given in (21 1, vt = ip' {^t4>t) — A^*) 

and hl^ is the t-th diagonal element of 



/2 
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Figure 1: Bias and RMSE for the estimators of a, /3 and 7 for = 0.98, constant precision 
model; ia (square), £p (circle), ire (triangle) and inaiveistsu). 
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Figure 2: Bias and RMSE for the estimators of a, /3 and 7 for kx = 0.84, constant precision 
model; ia (square), £p (circle), ire (triangle) and inaiveistsu). 
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Figure 3: Bias and RMSE for the estimators of a, j3, 7 and A for kx = 0.98, varying precision 
model; £a (square), ip (circle), £rc (triangle) and ^„aii>e(star). 
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Figure 4: Bias and RMSE for the estimators of a, (3, 7 and A for = 0.84, varying precision 
model; £a (square), ip (circle), £rc (triangle) and ^„oii>e(star). 
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(a) (b) 

Figure 5: Coverage of confidence intervals of a, /? and 7 for: column (a) kx = 0.98, constant 
precision model, and colimm (b) kx =0.84, constant precision model; (square), tp (circle), ire 
(triangle) and Inaivei^^)- 
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Figure 6: Coverage of confidence intervals of a, (3, 7 and A for: column (a) kx = 0.98, varying 
precision model, and column (b) = 0.84, varying precision model; £a (square), £p (circle), ire 
(triangle) and inaivei^^)- 



J with nit = (ptVt/lg'int)?, $ = diag{0 



1) 



in which M = diagjmi, . . . , m, 
(f)n} and W is an n X (pa + pp) matrix with the t-th row given by (z^, w^). Here, 
hat indicates that the unknown parameters are replaced by estimates. We suggest 
the use of the maximum pseudo-likelihood estimates since they performed well 
in our simulations and are computationally less demanding than the approximate 
maximum likelihood estimates. 

Plots of residuals versus observation indices are not always suitable for de- 



tecting lack of fit when measurement errors are present; see, for instance. Fuller 



(1987, Section 2.2.2), Carroll and Spiegelman (1992) and Buonaccorsi (2010 



Section 4.7). |Fuller| ( |T987l Section 2.2.2) suggests constructing plots of residuals 
versus consistent estimates of the expected value of Xt given wt- |Atkinson| ( |1985[ ), 
on the other hand, suggests the use of simulated envelopes in normal probability 
plots to facilitate their interpretation. The use of standardized weighted residuals 
plots proposed here will be illustrated in the next section. 



6. Real data application 

We now illustrate our results in the dataset described in Section [T] Our aim 
is to model RPP (relative pulse pressure) using BMI (body mass index) and 
LDL (low density lipoprotein cholesterol level) as covariates. The total choles- 
terol level, TC, is easier to be accessed than LDL, and provides a measure of 
LDL plus unknown quantities of other components as triglycerides and high den- 
sity lipoprotein. Here, we shall consider TC as a surrogate for LDL. 

We shall assume that yi,y2, ■ ■ ■ ,yn are independent observations of the PPR 
of the n individuals in the sample, and such that yt follows a beta distribution with 
mean fit and precision parameter (pt, with 

log[/it/(l - /it)] = ao + aiBMIt + /3LDLt, 

log(0t) = 70 + 7iBMIt + ALDLt, 

TQ = To + nLDU + et, (22) 

LDLt ~ N{f,,,al), et'^ N{0,a',), 

for t = 1, . . . , 182. Also, LDLt and e^, for t,t' = 1, . . . ,n are assumed to be 
independent Here, the interest parameter vector is 6 = (ao, ai, (3, 70, 71, A)^. 

For this particular dataset, the values of LDL and TC are available for all the 
individuals in the sample. We then use these data to mimic a situation where both 
the true and surrogate covariates are observed for some but not all individuals in 
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the sample. We randomly selected a subsample of 21 individuals for which we 
consider the corresponding observed values of LDL and TC; for the remaining 
individuals, only the observations on TC are considered in the analysis. A scatter 
plot of TC versus LDL for the selected individuals (not shown) suggests a clear 
approximate linear tendency. 

It is possible to estimate tq, ti and from (22), since we have observations 
on both TC and LDL for some individuals. We obtained the following estimates: 
To = 0.7351, Ti = 1.062 and a"^ = 0.030. To estimate the parameters of interest 
we used the approximate maximum likelihood, maximum pseudo-likelihood, re- 
gression calibration and naive methods. Also, as a gold standard to compare with 
these methods, we fitted a beta regression model in which LDL is used as the true 
covariate, measured without error, for all individuals. 

Table [T] shows the estimates, standard errors and p- values for the parameters 
of interest. All the approaches produce similar inferences on the parameters of 
the mean submodel (ao, cn and (3) and on the intercept (70) and the coefficient 
7i of BMI, which is the covariate measured without error, in the precision sub- 
model. Inference on A, the coefficient of LDL (the covariate measured with error) 
in the precision submodel, varies depending on the approach being used. The 
gold standard and the approximate maximum likelihood and pseudo-likelihood 
approaches indicate that the null hypothesis "Hq : A = should be rejected at the 
5% nominal level (p-value= 0.013, 0.025, 0.022, respectively), while T-Lq is not 
rejected when the regression calibration and the naive methods are employed (p- 
value= 0.064, 0.064). In other words, at the 5% nominal level, the approximate 
maximum likelihood and maximum pseudo-likelihood approaches agree with the 
gold standard in that they indicate that the precision varies with LDL, unlike the 
regression calibration and naive methods. Finally, it can be noticed that the re- 
sults for the approximate maximum likelihood and maximum pseudo-likelihood 
methods are very close. 

We now use the standardized weighted residual presented in Section[5]to inves- 
tigate the presence of outliers or any indication of lack of fit. Figure|7]shows resid- 
ual plots for the model fitted using the maximum pseudo-likelihood approach. 
Figure vj^a) shows the plot of the residuals against predicted values of LDL, LDL, 
and Figure |7];b) shows a normal probability plot with simulated envelope. There 
is no indication of outliers or any apparent pattern. This indicates that the errors- 
in-variables model considered here fits the data well. 
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Table 1: Estimates, standard errors and p- values 



ivieuiOQ 


Parameter 


Estimate 


Standard error 


z Stat 


p- value 








0.124 


—2 855 


004 






-0.009 


0.004 


-2.250 


0.023 


Gold standard 


B 


0.107 


0.050 


2.140 


0.034 


'Yn 
/u 


5.905 


0.764 


7.729 


0.000 




/I 


—0.022 


0.024 


—0.916 


0.356 




A 


-0.729 


0.292 


-2.497 


0.013 




rvr\ 


-0.366 


0.113 


—3.239 


0.006 






-0.009 


0.004 


—2.250 


0.024 




B 


0.118 


0.062 


1.903 


0.058 


Tn 


6.109 


0.908 


6.728 


0.000 




f 1 


-0.028 


0.028 


-1.000 


0.310 




A 


-0.751 


0.336 


-2.235 


0.025 






—0.366 


0.133 


—2.759 


0.006 




(X\ 


-0.009 


0.004 


-2.250 


0.013 


£p 


B 


0.118 


0.068 


1.735 


0.082 


Tn 
/u 


6.109 


0.994 


6.146 


0.000 




/I 


—0.028 


0.036 


—0.778 


0.434 




A 


-0.751 


0.329 


-2.283 


0.022 








n 1 ^4 

yj. 1 




yj.yjyjo 






-0.009 


0.004 


-2.250 


0.023 


f 




0.112 


0.062 


1.806 


0.068 


70 


5.945 


0.804 


7.394 


0.000 




71 


-0.028 


0.024 


-1.167 


0.239 




A 


-0.649 


0.351 


-1.849 


0.064 




ao 


-0.389 


0.144 


-2.701 


0.007 




ai 


-0.009 


0.004 


-2.250 


0.023 


^naive 


/3 


0.085 


0.047 


1.809 


0.068 


70 


6.136 


0.864 


7.102 


0.000 




71 


-0.028 


0.024 


-1.167 


0.239 




A 


-0.490 


0.265 


-1.849 


0.064 



7. Concluding remarks 

In this paper we proposed and studied errors-in- variables beta regression mod- 
els. We proposed three different estimation methods, namely, the approximate 

maximum likelihood, maximum pseudo-likelihood and regression calibration 
methods. We performed a Monte Carlo simulation study to compare the perfor- 
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(a) (b) 

Figure 7: Plot of standardized weighted residuals versus predicted values of LDL (a), and normal 
probability plot of standardized weighted residuals (b). 



mance of the estimators in terms of bias, root-mean- square errors and coverage 
of confidence intervals. Overall, we reached the following conclusions. First, ig- 
noring the measurement error may lead to severely biased inference. Second, the 
regression calibration approach is very simple and is reliable for estimating the 
parameters of the mean submodel when the measurement error variance is small. 
Third, the approximate maximum likelihood and maximum pseudo-likelihood ap- 
proaches perform well, the later being less computationally demanding than the 
former. We, therefore, recommend the maximum pseudo-Ukelihood estimation 
for practical applications. We also proposed a standardized weighted residual for 
diagnostic purposes. All our results were illustrated in the analysis of a real data 
set. 
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